home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / rvt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  7.7 KB  |  286 lines
  1. /*
  2.    ratio function:
  3.    sp = (usp1 - usp2) / (usp1 + usp2)
  4. */
  5.  
  6. #include <stdio.h>
  7. #include <spec.h>
  8. #define TMPFILE "global.def"
  9.  
  10. float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
  11. int   flg_v;
  12.  
  13. help()
  14. {
  15.   printf("ratio function: sp = (a*usp1 - b*usp2) / (a*usp1 + b*usp2)\n");
  16.   printf("rvt [-s1 file1] [-s2 file2] [-o file] [-ba1 n] [-ba2 n] [-last n]\n");
  17.   printf("options:\n");
  18.   printf("   -s1 file1    0 degree spctrum instead of USP1\n");
  19.   printf("   -s2 file2    90 degree spectrum instead of USP2\n");
  20.   printf("   -o  file     output spectrum instead of stdout\n");
  21.   printf("   -ba1 n       background for 0 degree spectrum\n");
  22.   printf("   -ba2 n       background for 90 degree spectrum\n");
  23.   printf("   -last n      takes arithmetic mean of last n channels as\n");
  24.   printf("                background. Default is, to find the background\n");
  25.   printf("                by fitting an exponential curve to the spectrum\n");
  26.   printf("   -v           verbose mode. Prints fitted background\n");
  27.   printf("   -def file    set spectra definition file other then global.def\n");
  28.   printf("\n(C) Rainer Kowallik\n");
  29.   exit(0);
  30. }
  31.  
  32. float sd(x,y)
  33. float x,y;
  34. {
  35. if(y!=0.0) return((float) (x / y));
  36. return(0.0);
  37. }
  38.  
  39. main(argc,argv)
  40. int argc;
  41. char *argv[];
  42. {
  43. int  n,i,max,flg_ba,lastn,ba1,ba2,xbeg,xend;
  44. char z[80],comment[80],nam1[80],nam2[80],defnam[80];
  45. float x,alpha,sum1,sum2;
  46. FILE *fp;
  47.  
  48.    strcpy(nam1,"USP1");
  49.    strcpy(nam2,"USP2");
  50.    strcpy(defnam,TMPFILE);
  51.    flg_ba=0;            /* defaults to exponential fit */
  52.    ba1 = -1; ba2 = -1;
  53.    flg_v=FALSE;
  54.  
  55.    if(checkopt(argc,argv,"-s1",z)) strcpy(nam1,z); 
  56.    if(checkopt(argc,argv,"-s2",z)) strcpy(nam2,z); 
  57.    if(checkopt(argc,argv,"-def",z)) strcpy(defnam,z);
  58.    if(checkopt(argc,argv,"-ba1",z)) {
  59.       ba1=atoi(z); flg_ba=1;
  60.       ba2=ba1;
  61.    }
  62.    if(checkopt(argc,argv,"-ba2",z)) {
  63.       ba2=atoi(z); flg_ba=1;
  64.       if(ba1==-1) ba1=ba2;
  65.    }
  66.    if(checkopt(argc,argv,"-last",z)) {
  67.       flg_ba=2; lastn = atoi(z);
  68.    }
  69.    if(checkopt(argc,argv,"-v",z)) flg_v=TRUE;
  70.  
  71.    spc=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  72.    err=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  73.    tim=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  74.    usp1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  75.    usp2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  76.    uer1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  77.    uer2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  78.  
  79. /* -----------------------------------------------
  80.    reading spectrum 
  81.    ----------------------------------------------- */
  82.  
  83.    max=readspec(nam1,usp1,uer1,tim,comment);
  84.    max=readspec(nam2,usp2,uer2,tim,comment);
  85.  
  86. /* -----------------------------------------------
  87.    generating error spectrum 
  88.    ----------------------------------------------- */
  89.    for(n=0;n<max;n++) {
  90.       uer1[n] = sqrt(usp1[n]);
  91.       uer2[n] = sqrt(usp2[n]);
  92.    }
  93.  
  94. /* -----------------------------------------------
  95.    reading spectra definition file
  96.    ----------------------------------------------- */
  97.  
  98.    fp=fopen(defnam,"r");
  99.    while(!feof(fp)) {
  100.       fscanf(fp,"%d %d %d %d\n",&i,&i,&xbeg,&xend);
  101.    }
  102.    if(flg_v) printf("xbeg = %d, xend = %d\n",xbeg,xend);
  103.    fclose(fp);
  104.  
  105.  
  106. /* --------------------------------------------
  107.    finding correct background
  108.    -------------------------------------------- */
  109.  
  110.    switch(flg_ba) {
  111.    case 0:
  112.       ba1=expbafit(usp1,uer1,xbeg,xend);
  113.       ba2=expbafit(usp2,uer2,xbeg,xend);
  114.       break;
  115.    case 1: /* allready done when parsing arguments */
  116.       break;
  117.    case 2:
  118.       ba1 = sumspc(usp1,max-lastn,max) / lastn;
  119.       ba2 = sumspc(usp2,max-lastn,max) / lastn;
  120.       break;
  121.    }
  122.    if(flg_v) printf("ba1=%d   ba2=%d\n",ba1,ba2);
  123.  
  124. /* --------------------------------------------
  125.    Subtracting background
  126.    -------------------------------------------- */
  127.    for(n=0;n<max;n++) {
  128.       usp1[n] = usp1[n] - ba1;
  129.       usp2[n] = usp2[n] - ba2;
  130.    }
  131.  
  132. /* --------------------------------------------
  133.    calculating normalization
  134.    -------------------------------------------- */
  135.    sum1 = (float) sumspc(usp1,xbeg,max);
  136.    sum2 = (float) sumspc(usp2,xbeg,max);
  137.    alpha = sd(sum1 , sum2);
  138.    if(flg_v) {
  139.       printf("sum usp1 = %f\nsum usp2 = %f\n",sum1,sum2);
  140.       printf("alpha    = %f\n",alpha);
  141.    }
  142. /* --------------------------------------------
  143.    Ratio function
  144.    -------------------------------------------- */
  145.    for(n=0;n<max;n++) {
  146.  
  147.       usp2[n] = alpha * usp2[n];
  148.       x = usp1[n] + usp2[n];
  149.       spc[n] = sd(usp1[n] - usp2[n] , x);
  150.  
  151.       sum1 = sd(2*usp2[n] , x*x);
  152.       sum2 = sd(2*alpha*usp1[n] , x*x);
  153.       sum1 = sum1 * sum1;
  154.       sum2 = sum2 * sum2;
  155.  
  156.       err[n] = sqrt((sum1 * usp1[n]) + (sum2 * usp2[n]));
  157.    }
  158. /* -------------------------------------------
  159.    now we can go and save the new spectra 
  160.    ------------------------------------------- */
  161.  
  162.    strcpy(z,comment);
  163.    n=instr("|",comment);
  164.    if(n>0) midstr(z,comment,0,n-1);
  165.    strcat(z," | RVT");
  166.    writespec("",spc,err,max,2,z);
  167.    free(uer1); free(uer2); free(usp1); free(usp2);
  168.    free(spc); free(err); free(tim);
  169.    exit(0);
  170. }
  171.  
  172. sumspc(spc,a,z)
  173. float spc[];
  174. int a,z;
  175. {
  176. int i,erg;
  177.  
  178.    erg=0;
  179.    for(i=a;i<z;i++) erg=erg+spc[i];
  180.    return(erg);
  181. }
  182.  
  183. /* ---------------------------------------------------------
  184.    Liner regression after logarithm of spectrum
  185.    This is only to calculate the chi**2, all
  186.    other values are lost.
  187.    UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
  188.    I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
  189.    --------------------------------------------------------- */
  190.  
  191. float linreg(spc,err,n,a,z)
  192. float spc[],err[];
  193. int n,a,z;
  194. {
  195. int i;
  196. float ysoll,x,y,sumx,sumxx,sumy,sumyy,sumxy,nges,
  197.       steigung,offset,erg,
  198.       *ylog;
  199.  
  200. ylog=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  201.  
  202. nges = (float) (z-a);    /* this we know for shure */
  203. sumx = 0.0;
  204. sumxx = 0.0;
  205. sumy = 0.0;
  206. sumyy = 0.0;
  207. sumxy = 0.0;
  208.  
  209. for(i=a;i<z;i++) {
  210.    x = (float) i;
  211.    y = spc[i];
  212.    y = y - ((float) n);
  213.    if(y<0.1) y=1.0;                 /* trap overflow on log */
  214.    y = (float) log((double) y);     /* linearization of exp function */
  215.    ylog[i] = y;                     /* store for later use */
  216.    sumx = sumx + x;
  217.    sumxx = sumxx + (x*x);
  218.    sumy = sumy + y;
  219.    sumyy = sumyy + (y*y);
  220.    sumxy = sumxy + (x*y);
  221. }
  222. sumx = sumx / nges;
  223. sumy = sumy / nges;
  224. sumxx = sumxx - (nges * sumx * sumx);
  225. sumxy = sumxy - (nges * sumx * sumy);
  226. sumyy = sumyy - (nges * sumy * sumy);
  227. steigung = sumxy / sumxx;
  228. offset = sumy - (steigung * sumx);
  229.  
  230. /* calculating chi**2 */
  231.  
  232. erg=0.0;
  233. for(i=a;i<z;i++) {
  234.    x = (float) i;
  235.    y = spc[i];
  236.    y = ylog[i];
  237.    ysoll = (steigung * x) + offset;
  238.    y = y - ysoll;
  239.    erg = erg + (y * y); 
  240. }
  241. free(ylog);
  242. return(erg);
  243. }
  244.  
  245.    
  246. /* ---------------------------------------------------------
  247.    Fitting background to exponential function.
  248.    This is done, by calculating the chi**2 with
  249.    linear regression for the logarithm of the spectrum.
  250.    We than look for a minimum of chi**2 for different
  251.    backgrounds.
  252.    UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
  253.    I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
  254.    --------------------------------------------------------- */
  255. expbafit(spc,err,a,z)
  256. float spc[],err[];
  257. int a,z;
  258. {
  259. int i,n,bamin,bamax,babest,inc,xraster,xdepth;
  260. float chisq,chibest;
  261.  
  262. xraster = 10;
  263. xdepth = 5;
  264. chibest = 1E10;
  265. bamax = sumspc(spc,z-10,z) / 10 ; /* get maximum background */
  266. bamin = bamax / 2;                /* get minimum background */
  267. for(i=1;i<xdepth;i++) {
  268.    inc = (bamax - bamin) / xraster;
  269.    if(inc<1) break;
  270.    n = bamin; while(n<=bamax) {
  271.       chisq = linreg(spc,err,n,a,z); /* just calculate chi**2 */
  272.       if(chisq < chibest) {
  273.      chibest = chisq;
  274.      babest = n;
  275.      if(flg_v) {
  276.      printf("new best: background = %d   chi**2 = %f\n",babest,chibest);
  277.      }
  278.       }
  279.       n = n + inc;
  280.    }
  281.    bamin = babest - inc;
  282.    bamax = babest + inc;
  283. }
  284. return(babest);
  285. }
  286.